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ABSTRACT 

With the large sample of young 7 -ray pulsars discovered by the Fermi Large Area Telescope (LAT), population synthesis has become 
a powerful tool for comparing their collective properties with model predictions. We synthesised a pulsar population based on a 
radio emission model and four 7 -ray gap models (Polar Cap, Slot Gap, Outer Gap, and One Pole Caustic). Applying 7 -ray and radio 
visibility criteria, we normalise the simulation to the number of detected radio pulsars by a select group of ten radio surveys. 

The luminosity and the wide beams from the outer gaps can easily account for the number of Fermi detections in 2 years of observa- 
tions. The wide slot-gap beam requires an increase by a factor of ~ 10 of the predicted luminosity to produce a reasonable number of 
7 -ray pulsars. Such large increases in the luminosity may be accommodated by implementing offset polar caps. The narrow polar-cap 
beams contribute at most only a handful of LAT pulsars. Using standard distributions in birth location and pulsar spin-down power 
(. E ), we skew the initial magnetic field and period distributions in a an attempt to account for the high E Fermi pulsars. While we 
compromise the agreement between simulated and detected distributions of radio pulsars, the simulations fail to reproduce the LAT 
findings: all models under-predict the number of LAT pulsars with high E, and they cannot explain the high probability of detecting 
both the radio and 7 -ray beams at high E. The beaming factor remains close to 1.0 over 4 decades in E evolution for the slot gap 
whereas it significantly decreases with increasing age for the outer gaps. The evolution of the enhanced slot-gap luminosity with E 
is compatible with the large dispersion of 7 -ray luminosity seen in the LAT data. The stronger evolution predicted for the outer gap, 
which is linked to the polar cap heating by the return current, is apparently not supported by the LAT data. The LAT sample of 7 -ray 
pulsars therefore provides a fresh perspective on the early evolution of the luminosity and beam width of the 7 -ray emission from 
young pulsars, calling for thin and more luminous gaps. 
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1. Introduction 

After the radio detection of the first pulsar signal in 1967 
(Hewish et al., 1968), a pulsar magnetosphere model was for- 
mulated by Goldreich & Julian (1969). A direct consequence of 
the Goldreich & Julian model is the establishment of a magne- 
tospheric charge density that creates a force-free pulsar magne- 
tosphere. However, such a magnetosphere has no electric field 
along the magnetic field to accelerate charges and produce 7 - 
rays. 

The detection, a few years later, of pulsed emission at y- 
ray energies from the Crab (McBreen et al., 1973) and Vela 
(Thompson et al., 1975) pulsars, and the detection of four more 
y-ray pulsars by Thompson et al. (1994) established that pul- 
sars accelerate particles to energies of at least a few TeV sug- 
gesting that there are magneto spheric regions where the charge 
density departs from that of Goldreich & Julian, locally vio- 
lating the force-free condition and allowing particle accelera- 
tion. These regions were identified in two magneto spheric zones. 
In the inner magnetosphere, acceleration can take place both 
above the polar cap and in the slot gap , which extends to high- 
altitude along the last open magnetic field lines. In the outer 
magnetosphere, the outer gap extends from the null charge sur- 


face to the light cylinder. These gap regions correspond to three 
models: the low-altitude slot-gap model, hereafter Polar Cap 
(PC, Muslimov & Harding (2003)), the Slot Gap model (SG, 
Muslimov & Harding (2004)), and the Outer Gap model (OG, 
Cheng et al. (2000)). 

In the polar-cap model the emission comes from a region 
close to the neutron star (NS) surface and well confined above 
the magnetic polar cap. Charged particles from the neutron star 
are initially accelerated in the strong electrostatic field gener- 
ated by a departure from the Goldreich- Julian charge density 
(Arons & S char lemann, 1979). Aided by inertial frame drag- 
ging (Muslimov & Tsygan, 1992), pulsars emit high energy pho- 
tons by curvature radiation (CR) and inverse Compton scattering 
(ICS). The most energetic of these photons reach threshold for 
electron-positron pair production in the strong magnetic field at 
a Pair Formation Front (PFF), above which the secondary pairs 
can screen the electric field in a short distance. The pairs, pro- 
duced in excited Landau states, emit synchrotron photons which 
trigger a pair cascade with high multiplicity. A small fraction of 
the pairs is actually accelerated. The pair plasma likely estab- 
lishes force-free conditions along the magnetic field lines above 
the PFF, as well as radiate 7 - rays. Over most of the polar cap, 
the PFF and 7 - ray emission occurs well within a few stellar radii 
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of the NS surface. The main contribution to the y- ray emission 
comes from CR from the pairs moving upward. Since the CR in- 
tensity scales with the magnetic field lines curvature, it decreases 
from the polar cap edge toward the magnetic axis, conferring to 
the emission beam the structure of an hollow cone. 

The slot-gap emission is generated from the same polar cap 
electromagnetic pair cascade near the boundary of the closed 
magnetic field lines region where the parallel electric field E\\ — > 
0 and the PFF rises to higher altitude. Here electrons are ac- 
celerated over longer distances to produce the pair cascade. A 
narrow gap, the slot gap , is formed along the closed magnetic 
field surface where the PFF is never established, and electrons 
continue to be accelerated and radiating y- rays by self-limited 
curvature radiation into the outer magnetosphere. The resulting 
hollow beam is much broader and less collimated near the mag- 
netic axis than the lower- altitude PC emission (see Section 7). 

The outer gaps are vacuum regions characterised by a strong 
electric field along the magnetic field lines (Holloway, 1973; 
Cheng et al., 1976) above the null charge surface. Two outer 
gap regions (Cheng et al., 1976; Romani & Yadigaroglu, 1995; 
Cheng et al., 2000; Hirotani, 2006) can exist in the angular 
velocity-magnetic momentum plane, one for each pole. In the 
physical OG model, in the case of a non-aligned rotator, the 
gap region closer to the pulsar surface is more active than the 
other gap further away from the surface due to the pair produc- 
tion screening operating more efficiently at lower altitude. In the 
OG model a charge-deficient region forms in the outer magne- 
tosphere above the null charge surface where a charge-separated 
flow is formed. The induced electric field accelerates pairs radi- 
ating y-ray s in a direction tangent to the B lines. The y - ray pho- 
tons interact with thermal X-rays from the NS surface to produce 
pairs on field lines interior to the last open field line. The pair 
formation surface screening the electric field defines the interior 
surface of the gap. 

More than 2000 pulsars are listed in the ATNF database 
(Manchester et al., 2005), most of which were first observed at 
radio wavelength. We employ the following ten selected pul- 
sar radio surveys in this study: Molonglo2 (Manchester et al., 
1978), Green Bank 2 & 3 (Dewey etal., 1985; Stokes et al., 
1985), Parkes 2 (70 cm) (Lyneetal., 1998), Arecibo 2 & 3 
(Stokes et al., 1986; Nice et al., 1995), Parkes 1 (Johnston et al., 
1992), Jodrell Bank 2 (Clifton & Lyne, 1986), Parkes Multi- 
beam (Manchester et al., 2001) and the extended Swinburne sur- 
veys (Edwards et al., 2001; Jacoby et al., 2009). For these, the 
survey parameters are known with a high accuracy and they 
cover the largest possible sky surface while minimising the over- 
lapping regions. 

The advent of the LAT telescope on the Fermi satellite 
(Atwood et al., 2009) led to a drastic increase in the number of 
y-ray pulsars. After three years of observations the LAT detected 
about 106 pulsars, more than doubling the number of detections 
listed in the first pulsar catalog (Abdo etal., 2010) leading to 
the discovery of two well defined y- ray pulsar populations con- 
sisting of 31 millisecond pulsars, and 75 young or middle aged 
isolated, normal pulsars. To study and compare the collective 
properties of the LAT normal isolated pulsars and investigate the 
emission mechanisms that best explain the observed emission, 
we synthesised a pulsar population incorporating four important 
high-energy radiation gap models. The simulation takes into ac- 
count the axisymmetric structure of our Galaxy and is designed 
to match the known characteristics of the group of older radio 
pulsar population than the younger group of pulsars sampled in 
y-rays. Four y-ray emission gap models have been assumed: the 
previously described Polar Cap (PC), Slot Gap (SG), and Outer 


Gap (OG), and a variation of the OG, hereafter the One Pole 
Caustic (OPC) (Romani & Watters, 2010; Watters et al., 2009) 
that differs from the OG in the energetics. We model the radio 
emission at two different frequencies, 1400 MHz and 400 MHz 
(Gonthier et al., 2004; Harding et al., 2007), comparing simu- 
lated radio fluxes with the flux thresholds of existing surveys. 

The outline of this paper is as follows. In Sections 2 and 
3, we describe the neutron star characteristics and evolution. In 
Sections 4, 5, and 6, we give a brief overview of the radio lu- 
minosity computation, y-ray gap widths, and y-ray luminosities 
computations. Sections 7 and 8 describe the pulsar light-curve 
and flux computation. Section 9 reviews the radio and y-ray 
pulsar visibility calculations. We present the results in the final 
Section 10. 


2. Neutron star characteristics 

The neutron star mass, radius, and moment of inertia used 
in this paper have been chosen according to the experimen- 
tal mass measurements in binary NS-NS systems, X-ray bi- 
naries, and NS-white dwarf binaries shown in Figure 3 of 
Lattimer & Prakash (2007). 

The assumed NS mass and radius are M N s = 1.5 M Q and 
7 ? n s = 13 km. The mass value lies between the weighted av- 
erage and average values of X-ray and white dwarf-NS bina- 
ries estimates and, with the 7? NS = 13 km, represent, a possible 
solution for the EOS that describe the NS interior (Figure 2 of 
Lattimer & Prakash (2007)). 

The moment of inertia of a NS is evaluated by Equation 
35 of Lattimer & Prakash (2007). For the 13 km radius and the 
1.5 M q mass of our standard NS, we obtain / ~ 1.8 x 10 38 kg m 2 . 
Because of the uncertainty on the mass and radius estimates, this 
value has an uncertainty of about 70%. 

For each simulated NS we have generated a value of the mag- 
netic obliquity a (angle between the pulsar rotation and mag- 
netic axes) and of the observer line of sight £ (angle between the 
pulsar rotation axis and the observer line of sight). After the su- 
pernova explosion that generates the neutron star, the magnetic 
axis a has equal probability to point in any direction of a 3 di- 
mensional space. This is also true for the observer line of sight 
direction £ with respect to the pulsar rotational axis. The a and 
£ distributions are isotropic. 

The spin-down power E is defined as the rate with which the 
pulsar loses rotational kinetic energy, as 

9 • 9 — dL'rot oq P P ~ 3 

E = 4 n 2 IPP- 1 = — -21 ~ 7.1 x 10 39 -— — W. (1) 
d t 1 s/s Is 

The latter equation is based on the NS structure assumptions 
through the moment of inertia /. Since mass and radius are cho- 
sen inside intervals of allowed values, the E estimate is affected 
by an uncertainty of at least a factor 3. 

The choice of mass, radius, and moment of inertia formula- 
tion yields a moment of inertia value that is 1.5 times higher 
than in the ATNF catalog. This helps to reduce the discrep- 
ancy found between the simulated and observed E distributions 
(Section 10.3), while remaining well within the range of pa- 
rameters allowed by the binary data and equations of state in 
Lattimer & Prakash (2007). The choice of different values for 
mass and radius would also impact the range of the P and P 
distributions of the evolved pulsar population (Section 10.2). 
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3. Neutron stars at birth and their evolution 

We synthesised ~ 2.7 x 10 8 NSs with mass, radius, and moment 
of inertia as described in Section 2, and assuming a constant 
birth rate over the last 1 Gyr. It yields 2.5 x 10 6 isolated ordinary 
pulsars to the left of the radio death line (see below). In order 
to match the observed radio pulsars P and P distributions, an 
exponential magnetic field decay with a time scale of 2.8 x 10 6 
yr has been assumed (Gonthier et al., 2004). The choice of such 
a short timescale decay is justified by the need to slow down 
the birth population enough to reproduce the characteristics of 
the observed radio sample. It provides a simple mathematical 
solution to a more physical model of the rotational evolution of 
the NS, yet to be developed. For our study, since we are dealing 
with young ordinary pulsars, this choice has been checked not to 
affect the obtained results. 

The radio death line we used is defined as 

log/* < a + b logR (2) 

It is composed by three different segments (Story et al., 
2007; Zhang et al., 2000), each one refers to a specific period 
interval characterised by the following a and b values 


lack of high E objects (Section 10.3). The Po birth distribution 
has been derived from Pq and Bq by using the equation 


B s = 


3 c 3 Ins 


PP 


1/2 


( 4 ) 


This formulation includes no dependence on the magnetic 
obliquity alpha, as proposed by Ruderman & Sutherland (1975) 
for the spinning down of magnetospheres carrying current flows. 
More recently, Spitkovsky (2006) numerically showed for force- 
free magnetospheres that the spin down of orthogonal rota- 
tors is twice that of aligned rotators. In the non-ideal case of 
a magnetosphere accelerating charges to produce pulsed emis- 
sion the impact of alpha on E is still under discussion, so 
we chose for this paper the alpha independent prescription of 
Ruderman & Sutherland (1975). Hereinafter, all luminosities are 
given as a function of E to judge how the uncertainty on the spin- 
down rate propagates. 


3.2. Birth location and velocity in the Galactic plane 


P< 15ms a = -19.00 6 = 0.814 

15ms < P < 300ms a = -17.60 b = 1.370 (3) 

P > 300ms a = -16.69 b = 2.590 

3. 1. Birth spinning and magnetic characteristics 

The distribution of period at birth, Po, plotted in the right panel 
of Figure 1, follows a single gaussian of width 50 ms, centred 
at 50 ms, and truncated at 0 to avoid negative periods. The same 
distribution was adopted by Watters & Romani (201 1) on the ba- 
sis of radio luminosity arguments but it differs from the choice 
of Takata et al. (2011) who selected the birth period randomly in 
the range 20 < Pq < 30 ms. 



Fig. 1. Left: The assumed surface magnetic field distribution at 
birth. Right: The assumed spin period distribution at birth. 


The magnetic field birth distribution Bo shown in the left 
panel of Figure 1 has been built as the sum of two gaussians 
in log 10 Bo [Tesla], both 0.4 in width, respectively centred at 
8.5 and 9.1, and with an amplitude ratio of 1:7/12. Our choice 
represents a compromise between that of Watters & Romani 
(2011), a single gaussian centred at 8.65 and width 0.3, and 
the Takata et al. (2011) one, a single Gaussian centred at 8.6 
and width 0.1. The high-i?o Gaussians provide energetic pulsars 
when evolved. 

Both the Po and Bo distributions have been optimized a pos- 
teriori to obtain, after evolution, a simulated pulsar sample as 
close as possible to the observed one by minimizing the observed 


To follow the dynamical evolution of the pulsars in the Galactic 
reference frame, we synthesised their birth position x, y, z in the 
Galaxy as well as their kick velocity and direction. 

We emulated the distribution of the NS progenitors by us- 
ing the location of the HII regions in the Galaxy. The latter are 
good tracers of massive stars because O-B stars are required to 
ionise the hydrogen bubbles. For the number density of pulsars 
at birth as a function of Galactocentric distance, we used the HII 
region profile recently obtained by Bania et al. (2010) from ra- 
dio observations that can probe HII regions to large distance with 
little absorption. Figure 2 shows the comparison between the 
Paczynski (1990) birth distribution used in earlier publications 
(Gonthier et al. (2004); Takata et al. (201 1)) and the HII region 
profile used here. Both distributions extend from the Galactic 
centre up to 40 kpc and have been normalised to 1 over the 
Galaxy. 



Fig. 2. Surface density of the new born neutron stars. The dashed 
curve represents the Paczynski distribution (Paczynski, 1990), 
while the adopted one following the distribution of radio HII 
regions, is shown as a solid curve. Both curves are normalised. 
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We assume that all the NSs are bom in the Galactic disk, 
with an exponential thin disk distribution with a scale height of 
50 pc (consistent with Watters & Romani (2011) that adopted an 
exponential thin disk with a 75 pc scale height) and with a sur- 
face density distribution defined in Figure 2. Due to the large 
supernova kick velocity, the neutron stars evolve quickly out 
of the plane of the Galaxy. The assumed kick velocity distribu- 
tion is the same as in Watters & Romani (2011) and Takata et al. 
(2011). It is described by a Maxwellian distribution, charac- 
terised by a mean of 400 km s -1 and a width of 256 km s -1 
(Hobbs et al., 2005). 

3.3. Evolution 

We have evolved both the pulsar position and velocity in 
the Galactic gravitational potential (described in Paczynski 
(1990) and Gonthier et al. (2002) Equations 17, 18, and 19, and 
Takata et al. (201 1)). The spin characteristics have been evolved 
to the present time assuming a magnetic dipole. 

The simulated pulsar population at birth is shown, in red, in 
the P-P diagram of Figure 3. Following Gonthier et al. (2002), 
by knowing the analytical expression for B{t) = /(i?o, 0? it is 
possible to follow the evolution of the spin parameters from the 
birth time to to the present time t p . The magnetic decay is de- 
scribed by 

B(t) = B 0 , 8 <T' /td (5) 

where td = 2.8 Myr is the decay timescale, and Bo ,8 is the birth 
magnetic field in units of 10 8 Tesla. 

Assuming magnetic dipole spin-down and initial period Po , 
the period and the period first time derivative at the present time 
can be obtained from Equations 7 & 8 of Gonthier et al. (2002). 
The simulated pulsar population after evolution is shown, in 
blue, in Figure 3. 



Period (s) 

Fig. 3. P-P diagram of the pulsar population at birth (in red), and 
the population evolved to the present (blue). 


4. Radio emission model 


After evolving the neutron stars in the Galactic frame, values 
of the radio dispersion measure (DM), and the radio scatter- 
ing measure (SM), are assigned to each star using the NE2001 
model (Cordes & Lazio, 2001). The sky temperature at 408 MHz 
(T sk y,408) for each star is obtained using the all-sky map from the 
study of Haslam et al. (1982). 

The empirical radio emission model we have implemented in 
our simulations follows the work of Gonthier et al. (2004) and 
Harding et al. (2007). We assume that the radio beam is com- 
posed of a core component originating relatively near the neu- 
tron star surface and a conical component radiated at higher alti- 
tude, both centered on the magnetic axis in the co-rotating frame. 
The adopted form of this model is similar to that proposed by 
Arzoumanian et al. (2002), based on the work of Rankin (1983) 
and Kijak & Gil (2003) and modified to include frequency de- 
pendence by Gonthier et al. (2004). The total flux at a given fre- 
quency from the two components seen at angle 6 to the magnetic 
field axis is 

S(0, v) = F mK e^ lp ^ + F cone e- (e ^^ (6) 


where 


Fi(v) 


-(1 + a t ) / v \“i +1 U 
v 1 50 MHz) Q/D 2 ' 


(7) 


The index i refers to the core or cone, oil is the spectral index 
of the total angle and frequency integrated flux for each compo- 
nent, Lj is the component luminosity, and D is the distance to the 
pulsar. The total solid angles of the Gaussian beams describing 
the core and cone components are 


f^core _ TPcore ($) 

Q CO ne = 2 7l il2 W e e (9) 


where the latter can be written as 


f^cone - ^cone(^GHz) °' 26 ( 10 ) 


where vghz is the frequency expressed in Giga-Hertz. The fac- 
tor 6o) CO ne represents the portion of Q con e that is independent of 
the frequency and used later in Equation (14). The width of the 
Gaussian describing the core beam is 

Poore = 1.5° (-) (11) 

where P is the pulsar period in seconds. The annulus and width 
of the cone beam are 


0 = (1 - 2.63 d w )p cone (12) 

W e — ^wPcone ( 13 ) 

where S w = 0.18 (Gonthier et al., 2006), and 

Pcone = 1 . 24 °r^(£) °' 5 ( 14 ) 

is the radius of the open field volume at the emission altitude 
derived by Kijak & Gil (2003), and 

/ p \ 0 - 07 / p V 0.3 

^ 4 o (kf^) y (v - r °- 26 (i5) 

r K g is in units of stellar radius. The ratio of the core-to-cone peak 
flux r is expressed as 


r 


v 

n\ — 
\n 


*core ^cone 


-0.26 


(16) 
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and requires Qf C ore-»cone -0.26 = 0.9, a core -a cone = -0.64 where 
vi = 1MHz. Gonthier et al. (2006), who carried out a study of 
20 pulsars having three peaks in their average-pulse profiles at 
frequencies 400, 600 and 1400 MHz, found a core-to-cone peak- 
flux ratio 

I 10 4 I ('^ V ' 3 y-°' 9 

J7 AU (is/ uHz’ 

_ r core _ I v ' 

Fmne in3.3Z£r 1 -®v-0.9 

( lu Us/ GHz' 

It is consistent with the ratio of Arzoumanian et al. (2002) at pe- 
riods above about 1 s, and predicts that pulsars with P < 0.05 s 
are cone dominated. Such a picture is supported by the study of 
Crawford et al. (2001) who measured the polarisation of a num- 
ber of pulsars younger than 100 kyr, finding that they possess a 
high degree of linear polarisation and very small circular polar- 
isation, typical of cone beams. The luminosities of the core and 
cone components are 


P < 0.7s 


P > 0.7s. 


(17) 


where 


/cone — 


/radio 

1 + ( 1 M))’ 


/core — 


/radio 

1 + r 0 ' 


(18) 


1 / 1 + cr C0re 
G \ 1 + Qf C one 


( tUcone 
f^core 


v 0 - 26 

K 1000 


-..t^cone Q^core 

K 50 


, o .! 
1 


(19) 


r\ is evaluated from Equations 16 and 17, or C ore = -1.96, Of CO ne = 
-1.32, viooo = 1000 MHz and V 50 = 50 MHz, and 


/ p \-i / p \ 035 

Z/ radio = 2.805 x 10 9 j I -j-y- j mJy kpc 2 MHz 

as modified from Arzoumanian et al. (2002). 


( 20 ) 


5. PC & SG: particle luminosity and gap width 

5. 1. Particle luminosity 


The slot gap region is defined between the last open magnetic 
field line, defined by the colatitude 6o ^ (QR/cf(l)) l/2 , and the 
magnetic field line with a colatitude value (1 -wsg) where wsg is 
the SG width expressed in units of the dimensionless colatitude 
of a PC magnetic field line, f = Q/Qq. 

It is possible to define the emission component from the PC 
pair cascades along the PFF that forms on the inside surface 
of the SG by assuming that mono-energetic radiation is emit- 
ted tangent to field lines (Muslimov & Harding, 2003). wsg is 
a function of pulsar period, P, and surface magnetic field, i? N s 
(Muslimov & Harding, 2003). The photons from the polar cap 
pair cascade are emitted in the region defined by 1 - wsg- The 
luminosity of the SG from each pole is 


tSG 


r-ln nQ { 

’ I <#PC 
Jo Jo 0 ( 


pip, 77 )® (f, rj)r 2 sin 0d6 ( 21 ) 


0o (1 -wsg) 


where p(£, 77 ) and ®(f, 77 ) are the primary charge density and 
potential as a function of the emission altitude 77 = r/R^s 
and of f, in units of NS radius, and 0 pc is the magnetic az- 
imuthal angle. Using the expressions for ® and for p from 
Muslimov & Harding (2003), the PC particle luminosity (from 
the low-altitude SG) is 


L? = EwlaO ~ ^)M1 - *)(1 - f ) cos 2 a + 


3 

x(l - Wsg + YqM'Ig)^ 2 ) 1 ) 


tm I /(l) , 
H( 1 ) TfW 


( 22 ) 


sin 2 a] 


where E = Q 4 i?^ s 7^ s /(6c 3 /(l) 2 ) is the spin-down power, k = 
0.15 l 3 s/R 3 6 , I 2,8 is the NS moment of inertia in unit of 10 38 kg 
m 2 , Re is the NS radius R^s in unit of 10 6 m, 7 /( 77 ) is a relativis- 
tic correction factor of order 1 , firj) is the correction factor for 
the dipole component of the magnetic field in a Schwarzschild 
metric, and a is the pulsar obliquity (Muslimov & Tsygan, 1992; 
Harding & Muslimov, 1998). 

Using the equations for ® and for p from 
Muslimov & Harding (2004), the high-altitude SG particle 
luminosity from each pole can also be determined from 
Equation (21) as 

Ll G = £w 2 G //(l - + 1 (l - w SG + Tw 2 h )s (23) 

where = (1 - 377 / 477 i c ) 1/2 and r]\ c = nJR m = c/£ 1R NS - The 
parameters and !B , are defined as: 




+ 1 -/* 


1 \ 2 

+ cos a 

77lc/ 


<B=--H(l)H(ri) 


H(n c ) 

1 7(1) R 

/(l) 

77(1) 1 

y ’“/(>;.) -"] 

7] f(r])_ 


1/2 


% 


sin 2 a 


where tj lc = Rlc/^ns, fi = (1 - 0 . 75 t 7 /t 7 LC ), and tj c = 1.3. 
According to Muslimov & Harding (2004), the energies of the 
primary electrons in the SG quickly become radiation-reaction 
limited, with the rate of acceleration balancing the curvature ra- 
diation loss rate, resulting in 100% efficiency with L y = L S Q G in 
this case. 


5.2. Gap width 

In the SG model, the width of the slot gap wsg can be esti- 
mated as the magnetic colatitude where the variation in height 
of the curvature radiation PFF zo (in units of stellar radius) 
becomes comparable to a fraction A of the stellar radius R^s 
(Muslimov & Harding, 2003): 



In Equation 24, zo represents the dimensionless altitude, above 
the polar cap, of pair formation due to curvature radiation 


zq = 7 x 10 


-2 


d7/4 

0.1 


s 8 4 4 ^ 1/2 (i 


. £2)3/4 


(25) 


where P 0.1 = P/0.1 s, and is the magnetic field in units of 
10 8 Tesla. By solving numerically Equation 24 with zo defined 
in Equation 25, one obtains fsG for a specific pulsar. The wsg 
gap width value is then obtained as 


wsg = 1 - fsG- 


(26) 


The A parameter constrains both the energetics and emission 
pattern of the SG emission and impacts both the SG and PC 
luminosity (Section 10.6) and light-curve sharpness and shape. 
For large A values the light-curve peaks appear too sharp com- 
pared with the observed LAT profiles, therefore the slot gap is 
too narrow and not energetic enough to explain the observed 
LAT fluxes. On the other hand, smaller A values imply wider 
slot gaps, sufficiently luminous when compared with the obser- 
vations, but light-curve peaks too broad when compared with the 
observed ones. 
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As a result we compromise between the narrow light-curve 
structures and the y-ray luminosity through a reasonable radia- 
tion efficiency e y . We tried two different approaches to constrain 
A: one based on energetic arguments, and one based on the op- 
timisation of the expected light-curves for some of the LAT pul- 
sars. 

Since L y scales as w 3 SG x E and since we want the luminosity 
to be close to L y oc E l/2 (Abdo et al. (2010), First pulsar catalog) 
we need to have wsg E~ l/ 6 to obtain a reasonable agreement 

with the LAT data. The luminosity remains close to E 0 - 5 for all 
the tested A values, but favours A < 0.4 to explain the bright LAT 
pulsars. A good compromise is found for A = 0.35. 

One can calculate numerically the pair formation front shape 
for the P and B values of some of the best known pulsars, Crab, 
Vela, CTA1, and Geminga, to obtain an approximate wsg value 
(Muslimov & Harding (2003)). The results yield wsg, Crab =0.03, 
wsg, V ela = 0. 1, wsg,ctai=0.16, wsg, G eminga =0.3 for A values be- 
tween 0.02-0.6. 

In order to investigate how the pulsar light-curve changes as 
a function of A , we performed a fit to some LAT light-curves 
with the SG phase-plots (see Section 7), evaluated for a set of 
wsg values obtained for different A values. We studied the be- 
haviour of the best-fit likelihood value as a function of A for Vela, 
Crab, J1028-5820, J1048-5832, J2021+3651, and J2229+61 14. 
For all the studied pulsars, in the A < 0.4 range that allows bright 
enough pulsars, the maximum-likelihood value presents a local 
maximum between 0.2 and 0.4. This result is consistent with the 
A estimate obtained from the luminosity study and the pair for- 
mation front evaluation from Crab, Vela, CTA1, and Geminga. 

In this paper, we set T=0.35. This value reproduces the bulk 
of the light-curve structure of the observed objects and yields a 
reasonable estimate of the SG luminosity. In choosing A, we put 
more emphasis on matching the sharply peaked light-curves of- 
ten recorded by the LAT than on achieving bright luminosities. 
This selection of A was driven by the need to preserve realis- 
tic beam patterns (thus their brightness and visibility across the 
beam) and is a key assumption that contains the results of our 
population studies. We mitigated the low SG gamma-ray lumi- 
nosities by using a radiative efficiency greater than 1 as discussed 
in section 8.2. 


6. OG and OPC: particle luminosity and gap width 

6.1. Gap width 

To determine the gap width, we consider two different prescrip- 
tions. The first one (Watters et al., 2009) simply assumes that the 
gap width is equal to the y-ray radiation efficiency. Because of 
the L y oc E 0 - 5 relation observed in the first LAT pulsar catalog 
(Abdo et al. (2010)) , the gap width should follow as 


Wo pc - 


/10 26 W 
l E 


1/2 


(27) 


Our second prescription follows the calculations of the self 
sustaining OG model presented in Zhang et al. (2004). In this 
formulation, the X-rays that trigger the pair production come 
from the bombardment of the NS surface by the full return cur- 
rent from the OG. The bright X-ray luminosity allows active 
OGs and y-ray emission for many old pulsars. The outer gap 
width across magnetic field lines is determined by computing the 
location of the pair formation surface. From Kapoor & Shukre 
(1998), the polar angle 6 C corresponding to the magnetic field 


line tangent to the light cylinder is: 


tan# c 


3 

4 tan a 


1 + (1 + ^ tan 2 a) 0 ' 5 


(28) 


with the light cylinder radius given by 


*L = 


r c 

sin 6 C ’ 


(29) 


Here r c is the distance between the pulsar and the point where the 
light cylinder is tangent to the magnetic field line corresponding 
to 0 C . The lower boundary of the outer gap is estimated from 
the null-charge surface, il • B = 0, that in two dimensions is 
described by (r m , 6- m ). By definition, the polar angle at the inner 
edge of the outer gap is 

tan0i n = ^3 tan a + V9 tan 2 a + 8 j . (30) 


The computation of r m is obtained from the relation 

sin 2 (0 - a) _ sin 2 (# c - a) 
r r c 

which results in 

R l sin 2 (6» in - a) 

n n = rr: 

sin0 c sin (0 C - a) 


(31) 


(32) 


The relation that defines the fractional OG size in this case is: 


woo = 5.2B- 4/7 P 26 ' 21 R- lon G((r), a) = /«r>, a). (33) 

where G(r, a) is a factor that is numerically solved for each pul- 
sar by taking into account the average distance (r) at which pri- 
mary y-rays are produced and along which magnetic field line 
they pair produce when they interact with an X-ray coming radi- 
ally from the NS surface. The average distance (r) is defined in 
Zhang et al. (2004) as 


a)rdr 
£" f((r), a)dr 


(34) 


where r max = min(r c , rb) and is the radius at which the frac- 
tional size of the outer gap stops to grow: /(rb, a)=l. 

A full calculation of the width of the OG radiating layer is 
complicated (Hirotani, 2006, 2008) since both the screening and 
the radiation occur in the same location. For this paper, we as- 
sume that this is an infinitely thin layer on the gap inner edge and 
that it is uniform in azimuth around the magnetic axis whereas 
Hirotani (2006) finds a significant azimuthal dependence. 


6.2. Particle luminosities 

The assumed gap width wopc defined in section 6.1 is not based 
on any physical prescription and is very different from the usual 
dependence luminosity oc(gap width) 3 (both SG and OG) based 
on the electrodynamics. 

The gap width luminosity is evaluated as 

L y , opc = wopc^sd- (35) 

In the OG case, from Zhang et al. (2004) and previous papers 
dealing with OG gap geometry, the total y-ray luminosity is 

L y , oc = w 3 0G ({r))E si (36) 

where wog is the fractional width of the gap at the average gap 
radius (r). 
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7. Phase-plot and light-curve generation 

7. 1. Assumptions and photon distributions 


To provide the y-ray emission pattern for each emission mech- 
anism, we used the geometric emission model from Dyks et al. 
(2004) based on the following assumptions: (i) the pulsar mag- 
netic field is dipolar and swept-back by the pulsar rotation (re- 
tarded potentials) (Deutsch, 1955), (ii) the y-ray emission is tan- 
gent to the magnetic field line and oriented in the direction of the 
accelerated electron velocity in the star frame. Relativistic aber- 
ration and time of flight delays are taken into account. 

In the computation of the emission pattern, the first step con- 
sists in localising the position of the magnetic field line from 
which the radiation is produced. Each field line is then divided 
into segments and for each segment the tangent direction and 
height with respect to the NS surface is evaluated. Since the 
emission gap is located, for each model, in a different magne- 
tospheric region, the emission patterns are obtained by selecting 
the segments corresponding to the gap position in each model. 
The y-ray emission is assumed to be uniform along the field lines 
in the co-rotating frame. The phase 0 of the pulsar emission is 
defined by the direction of the emitted photons with respect to 
the corotating frame. The result of this computation is the two- 
dimensional emission pattern in the plane ( 0 , 0 ), shown for each 
implemented emission model in Figure 4, which we refer to as 
a phase-plot. Figure 4 also shows the evolution of the emission 
pattern as a pulsar ages. 

To incorporate the radio emission geometry we modulate the 
field lines with the flux S ( 6 , v) given by Equation 6 . The differ- 
ential flux radiated from a bundle of field lines centred at open- 
volume coordinates (r, /) (see Dyks et al. (2004)) is 

dS ( 6 , v) = S i(6, v) sin 6 <7 0VC 6o r max — dv (37) 

N\ 

where Ni is the number of azimuthal divisions of each polar cap 
ring, r min = 0.1 and r max = 1 .0 are the lower and upper boundary 
of the emission region, and <7 0VC (= 0.1 for the radio phase-plot) 
is the spacing of the rings on the PC in open volume coordinates. 
For the SG model d OYC is adjusted to have 20 rings within the gap. 
The flux is assumed to be emitted at an altitude of 1 . 8 R N s for the 
core component and at an altitude given by Equation 15 for the 
cone component. 

In the PC model, the emission profile in colatitude is in- 
finitely thin along the inner edge of the slot gap (with wsg de- 
fined in section 5.2) and the intensity of the emission along the 
field line, Tpc, exponentially decreases from the polar cap edges 
to the magnetic pole 


| ex p(T?)’ 

Tpc - \ 

l ex p(^f)’ * > 


(38) 


Here Sf = 2.5, <Xi n = 1.0, cr ou t = 2.0, and 5 is the curvilinear 
distance along the field line starting from the NS surface. Both 5 
and Sf are in unit of Rns- 

To model the emission component from primary electrons in 
the SG model, we assume that radiation is emitted along the field 
lines in the slot gap, up to altitude 77 = rj mSLX (where 77 = r/R NS ). 
We assume an emissivity distribution across the SG as: 


m») = (i -£) 

(39) 

where 0* = 0 at the center of the SG and 


WSG 

(40) 


Such a distribution, that peaks in the centre of the gap and de- 
creases to zero at the gap edges, follows from the 0 * distribution 
of the SG potential (Muslimov & Harding, 2004). 

For the OG/OPC model, we describe the emitting region as 
an infinitely thin layer along the inner surface of the gap. The 
radio and PC phase-plots show the hollow cone patterns centered 
on the magnetic pole, while the SG and OG phaseplots show the 
caustic emission patterns characteristic of outer magnetosphere 
emission (Dyks et al., 2004). 

7.2. Light-curve generation 

The {ff) phase-plot space has been sampled in 180 x 90 bins. 
Each bin t?(0, 0) of the phase-plot gives the number of photons 



Fig. 4. The top to bottom panels illustrate the y-ray emission pat- 
tern obtained for a young (left) and old (right) pulsar, respec- 
tively for the PC, SG, OG/OPC, and radio (core plus cone) mod- 
els. For the PC and Radio models, the time evolution is obtained 
for B - 10 8 T and period increasing from 30 to 1000 ms. In the 
SG and OG/OPC models, the time evolution is obtained for a 
gap width wsg increasing from 0.04 to 0.5 and wog/opc increas- 
ing from 0.01 to 0.4 respectively. All the plots are given for an 
obliquity a = 45°. 


per solid angle per primary particle that can be observed in the 0 
direction at the rotational phase 0 : 


77(0, 0) 


d5i ph 

sin 0(70(70 ’ 


(41) 


Each phase-plot is obtained for a specific set of pulsar parame- 
ters that define its magneto spheric structure: the spin period P, 
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the surface magnetic field B N s, and the magnetic obliquity a. For 
the studied models, the phase-plot has the following dependen- 
cies: 

«y, PC/Radio = f(P, B , 

« r ,SG = /(wsg,») and w S g = f(P,B) 
og/opc = /(w,or) and w 0 g/opc = f(P*B). 

For each emission model, we have evaluated phase-plots for 
a values, from 5° to 90°, with a step of 5°. For each a value, the 
phase-plots have been evaluated for 2 magnetic field values and 
9 spin period values for the PC and radio models, and for 16 gap 
width values in the SG and OG/OPC cases. The complete set of 
sampled parameters is listed in Table 1 . 

To obtain the light-curve of a given NS, with a particular 
set of P, B N s, a, and gap width parameters, we interpolated 
the phase-plots noted in Table 1. When comparing phase pro- 
files for a different set of parameters, typically one profile will 
be narrower than another one making it nontrivial to interpolate 
between them. We adopted a non-linear interpolation which ex- 
pands the narrower light-curve covering the smallest phase range 
up to the phase extent of the wider profile, then applies a linear 
interpolation, and contracting the expanded profile back down to 
the extent of the original parent profiles. This strategy preserves 
the thin peaks and high degree of modulation that characterises 
the pulsar emission profiles at radio and y-ray wavebands. 

8. Flux calculations 

8. 1. Phase-plot normalisation and energy flux 

Let us define L po i e as the radiative luminosity from each pole, 
either in the y-ray s or in the radio. Assuming a value for the 
primary particle production rate N e and the energy E of each 
photon, one obtains a radiation luminosity per phase-plot bin: 


8.2. Computations: gap width and energy flux 
We calculated the y-ray and radio light-curve for each pulsar of 

r-2n 

the sample, storing the value of the integral J Q «(£,bs> 0)<# f° r 
the flux computation. The width of the emission gaps is com- 
puted using Equation 26 for the SG, and equations 33 and 27 for 
the OG and OPC. Because the PC and SG models do not apply 
when the gap becomes too large (pair-starved gaps should then 
be used), the flux for gap widths larger than 0.5 has been set to 
0. Because no emission remains visible from the thin inner edge 
of OG/OPC gaps when the gap width exceeds 0.7 the flux for 
gap widths larger than 0.7 has been set to 0. So all the pulsars 
with a gap width above these threshold levels are assumed to not 
produce any y-ray emission. 



3 

RL RQ 


4 5 6 7 3 

log(char. age) (yr) 


4 5 6 7 

log(char. age) (yr) 


dL = N e En((p, if) sin £d£d(p = An((p , if) sin £d£d(p (42) 


where A is a proportionality constant. One can normalise the 
phase-plot to the total radiation luminosity over the two poles 
according to: 

J f*7T r*2jx 

sin £d£ I n(<p,£)d<p (43) 

0 Jo 

We define the specific intensity / as 


1 = 


dL_ 

la 


0 


An((p, f) sin i^di^dcp 
sin ^d^dcj) 


Mt,0- (44) 


It is now possible to obtain the average energy flux observed 
by an Earth observer for a line of sight 


(vFy) = 


Jq f(^obs>» (p)d(j) 

2^D 2 


(45) 


Here, D is the pulsar distance. 

From the Equations 44 and 45, we can write the average en- 
ergy flux observed at the Earth as: 


Lpole fo ” n ^ obs ’ $ >4 ' P 

nL>2 Jo T s * n f n(</>, £)d<f> 


(46) 


The latter Equation establishes the relation between the lumi- 
nosity derived in the framework of a given model, L po i e , and the 

rQ.71 

integral of the pulsar light-curve J Q n(£ 0 b s , 0)<i0, obtained, from 
the phase-plot, for f This is related to the beaming factor 

fa discussed in Section 10.5 


Fig. 5. Number density of the visible y-ray pulsars obtained for 
each model as a function of characteristic age and energy flux 
times the square of the pulsar distance. These parameters can all 
be measured from the observations. The linear grey scale satu- 
rates at 1 star/bin for the polar cap and 2.5 star/bin for the other 
models. The pink contours outline the region where simulated 
radio-loud y-ray pulsars are found (at 20% of the maximum den- 
sity). The pink and green lines show the data for the radio-loud 
and radio-quiet LAT pulsars, respectively. 


For the radio luminosity computation, L po i e , we have used 
Equation 20 to evaluate the total radio luminosity and Equation 
18 to evaluate the luminosities of each core and cone compo- 
nent. The y-ray luminosity has been obtained by scaling the 
particle luminosities derived in Equations 22, 23, 36, and 35, 
for the PC, SG, OG and OPC models respectively by using a 
radiative eficiency s y . The latter has been chosen to provide 
a good agreement between the observed and simulated S y D 2 
distributions as a function of characteristic age (S y is the pho- 
ton flux and D the pulsar distance). This distribution involves 
only readily observable quantities. The solution adopted for each 
model is shown in Figure 5. The choice of radiative efficien- 
cies are: epc = 1.0, 6 sg = 12.0, 6og = 1.0, and 6 opc = 0.5. 
The high value of 6 sg needed for the SG requires either a su- 
per Goldreich-Julian current or a stronger value of the accel- 
erating electric field in the gap compared to the original calcu- 
lation by Muslimov & Harding (2004). This is quite possible if 
the polar cap is slightly offset, i.e., non- symmetrical around the 


Pierbattista et al. 2012: Gamma-ray pulsar population 



B 

Tesla 

p 

milliseconds 

a 

Degrees 

Gap Width values 

PC 

10 s , 10* 

30, 40, 50, 75, 100 
300, 500, 750, 1000 

5-90 
5° step 

0.04, 0.06, 0.08, 0.1, 0.13, 0.16, 0.2, 0.225 
0.25, 0.275, 0.3, 0.34, 0.38, 0.42, 0.46, 0.50 

SG 

none 

none 

5-90 
5° step 

0.04, 0.06, 0.08, 0.1, 0.13, 0.16, 0.2, 0.225 
0.25, 0.275, 0.3, 0.34, 0.38, 0.42, 0.46, 0.50 

OG/OPC 

none 

none 

5-90 
5° step 

0.01, 0.025, 0.04, 0.05, 0.067, 0.084, 0.1, 0.2 
0.3, 0.4, 0.5, 0.53, 0.56, 0.59, 0.62, 0.65 

Radio 

10 8 , 10* 

30, 40, 50, 75, 100 
300, 500, 750, 1000 

5-90 
5° step 

none 


Table 1. Magnetic field, period, and gap width values for which the phase-plots have been evaluated for each emission model. The 
SG and OG/OPC emission patterns do not depend directly on the pulsar period and magnetic field. 


magnetic axis, as one expects from the shape of the magnetic 
field lines distorted by the stellar rotation. Harding & Muslimov 
(2011) show that this distortion leads to a larger pair multiplicity 
as well as an increased electrical field along the field lines, thus 
an enhanced luminosity. Offset polar caps can sustain the modest 
increase in particle energy that is required in the present popu- 
lation study to account for the flux and pulsar counts observed 
by the LAT without invoking a radiation efficiency larger than 
one. The offset polar cap prediction was not available at the time 
of the population synthesis work, so we keep here the original 
polar cap luminosity and 6sg = 1200%. 


8.3. Gamma-ray energy and photon flux 

In order to evaluate the photon flux from the energy flux and to 
compare it with the LAT sensitivity in photon flux, we need to 
assume an emission spectrum for the pulsars. The typical photon 
spectrum of a LAT pulsar is well fitted by a power-law with an 
exponential cutoff, like 


dN y 

~dE 


= k 



(47) 


In the first LAT pulsar catalog (Abdo et al., 2010), the distri- 
bution of the parameters T and E cvX can be described by two 
Gaussians: 


Gaussian X : mean = 1.97; variance = 0.18 
Gaussian Y : mean = 3.06; variance = 0.37 


(48) 


The spectral index T and the log^^cut) are defined as: 

6 = 0.5982 [rad] 

T = X cos 6 - Y sin 6 (49) 

log 10 (is cu t) = Xsin 6 + Y cos 6. 

The Gaussian widths, centroids, and correlation angle had been 
derived from the analysis of the spectral parameters measured 
for the 1st LAT pulsar catalogue. We took here the same values. 

The photon flux computation has been done using Equation 
46, 47, 48, and 49 by assuming that the luminosity is mainly pro- 
duced at photon energies >100 MeV. The choice of this thresh- 
old and the choice of radiative efficiencies to match the data in 
Figure 5 are linked. 


9. Gamma-ray and radio visibilities 

9. 1. y-ray pulsar visibility 

In order to select the simulated pulsars that could be detected by 
the LAT during two years of observation, we made use of the 
6 month pulsar visibility map published for the 1 st LAT pulsar 
catalogue (Abdo et al., 2010) and of the 1 year pulsar visibility 
of blind pulsar searches (Dormody et al., 201 1). The two maps 
have been used to estimate the y-ray detectability of the radio- 
loud pulsars (corresponding to the LAT radio selected objects) 
and the radio-quiet ones (corresponding to the LAT blind search 
objects) respectively. The maps give the minimum visible pho- 
ton flux S m i n .ph and have been obtained taking into account the 
real LAT observation time in the sky, the photon energy, and the 
effective collection area corrected for the different incidence di- 
rections. Since the sky survey mode for LAT observations has 
been continued after 6 months, the maps have been scaled to 2 
years as the square root of time for the radio-selected sensitivity 
map and linearly with time for the blind search sensitivity map. 
Very few pointed observations were programmed that would sig- 
nificantly alter the shape of the visibility map. Photons collected 
in survey mode largely dominate and the flux threshold for de- 
tectability is primarily limited by the intensity of the interstellar 
background. 

9.2. Radio visibility 

The synthesis of the population is not based on any assumed NS 
birth rate; we assume a flat star formation rate over the last 1 Gyr. 
Instead, the simulated sample was scaled to the real number of 
pulsars detected in the Galaxy. The scaling factor has been eval- 
uated by selecting all the ATNF radio pulsars present in a select 
group of ten surveys and comparing this number with simulated 
radio pulsars visible in the same region. We have generated a 
large enough population to reduce the Poisson fluctuations and 
to improve the statistics in the analysis results. 

We selected the simulated pulsars within the visibility cri- 
teria of 10 radio surveys from the ATNF database 1 for which 
the survey parameters are well known and that cover the largest 
possible sky surface while minimising the overlapping regions. 
These surveys are: Molonglo2 (Manchester et al., 1978), Green 
Bank 2 & 3 (Dewey et al., 1985; Stokes etal., 1985), Parkes 2 
(70 cm) (Fyne et al., 1998), Arecibo 2 & 3 (Stokes et al., 1986; 
Nice et al., 1995), Parkes 1 (Johnston et al., 1992), Jodrell Bank 
2 (Clifton & Fyne, 1986), Parkes Multi-beam (Manchester et al., 
2001) and the extended Swinburne surveys (Edwards et al., 


1 http ://www. atnf. csiro . au/research/pulsar/psrcat/ 
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^(KJy-l) 

<T S/N 

Tree (K) 

^surv (MHz) 

T’int (S) 

4 am p (ms) 

Av b (MHz) 

Av ch (MHz) 

Molonglo 2 

5.100 

5.4 

225 

408 

40.96 

40.0 

4.0 

4.0 

Green Bank 2 

0.886 

7.5 

30.0 

390 

136 

33.5 

16.0 

2.0 

Green Bank 3 

0.950 

8.0 

30.0 

390 

132 

2.0 

8.0 

0.25 

Parkes 2 

0.430 

8.0 

50.0 

436 

157.3 

0.6 

32.0 

0.125 

Arecibo 2 

10.9 1 1 

8.0 

100 z 

430 

39.3 

0.4 

0.96 

0.06 

Arecibo 3 

13.35 3 

8.5 

70.0 4 

430 

67.7 

0.5 

10.0 

0.078 

Parkes 1 

0.256 

8.0 

45.0 

1520 

157.3 

2.4 

320.0 

5.0 

Jodrell Bank 2 

0.400 

6.0 

40.0 

1400 

524.0 

4.0 

40.0 

5.0 

Parkes MB 

0.460 

8.0 

21.0 

1374 

2100.0 

0.25 

288.0 

3.0 

Swinburne 

0.427 

10.0 

21.0 

1374 

265.0 

0.25 

288.0 

3.0 


1 

2 
3 


4 


Computed using Gf- = 19 ~ < ° 1 ^ x 7 l s 19 ~ <i|) . 
Computed using T rec = 90 + 2.083 x 1 19 - 6\) 


Computed using Sf = 

Computed using T rec = 65 + 2.083 x 1 19 - d|). 


Table 2. Instrumental parameters of the radio surveys. For the Arecibo surveys we chose to adopt a more accurate definition for the gain and the 
receiver temperature that are function of the declination 6. Respectively from the left to the right column, are indicated: telescope gain divided by a 
system losses factor, minimum signal to noise detected, receiver temperature, central observation frequency, integration time, sampling time, total 
bandwidth, and channel bandwidth. 


2001; Jacoby et al., 2009). The ratio between the number of sim- 
ulated pulsars meeting the surveys’ visibility criteria and the 
number of objects actually detected is 


Sf 


^obs,10 surv 
^sim,10 surv 


(50) 


This is the factor we used to scale the simulated pulsar sample. 


9.2.1. Radio pulsar selection 

During a radio survey, the edges of the survey region are defined 
by the position of the radio-telescope beam centre. Nevertheless, 
because of the solid angle extension and complexity of the beam, 
it is possible to observe a pulsar slightly out of the declared sur- 
vey region. Thus, to say that all the pulsars observed during a 
survey fall inside the declared survey coordinates edges is not 
totally correct. The first parameter we re-evaluated for each sur- 
vey is the number of pulsars seen inside a given region. 

The second important parameter is the survey efficiency 6 surv . 
It is defined as a filling factor, e.g. the ratio between the actual 
solid angle covered by the radio telescope beam during the ob- 
servations, and the area within the declared survey boundaries. 
The survey efficiency can be considered as the probability of ob- 
serving a pulsar present in the survey region only if the parent 
spatial distribution is uniform. To evaluate the boundaries of the 
survey region and to define the survey efficiency we decided: 

1 . to slightly extend the sky survey boundaries in order to in- 
clude the largest number of pulsars actually detected by a 
survey, without changing too much the original boundaries 

2. to evaluate the detection flux threshold for each pulsar within 
a survey by scaling the Dewey formula (Dewey et al., 1985) 
with a free parameter, 6 DeW ey ? to match the observations. 

S min = ^Dewey X S threshold (51) 

where the threshold flux S threshold * s expressed by the Dewey 
formula 


c <TS/N[^rec + ^sky(A ^)] j W 

threshold = G -yjN^Bt V T^W' 


(52) 


The Dewey formula, or radiometer formula, takes into account 
the characteristics of a given radio telescope and detector as well 
as a pulsar period and direction to give the minimum flux the sur- 
vey would be able to detect. In Equation 52 ct S /n is the minimum 
signal to noise ratio taken into account, T rec is the receiver tem- 
perature, T s k y is the sky temperature at 408 MHz, G = Gain IP 
is the ratio between the radio telescope gain and the dimension- 
less factor p that accounts for system losses, N p is the number 
of measured polarisations, B is the total receiver bandwidth, t is 
the integration time, and P is the pulsar period. W is the effective 
pulse broadening, defined as 

W 2 = W 2 + r^ amp + t 2 dm + T^ cat + r 3 ailDM . (53) 

Here, Wo is the intrinsic pulse width (Duty Cycle), r samp is a low- 
pass filter time constant applied before sampling (when this pa- 
rameter is unknown, a value equal to twice the sampling time has 
been used), r D m is the pulse smearing due to the DM over one 
frequency interval Av, and T sca t is the pulse broadening due to in- 
terstellar scattering (Dewey et al., 1985). The dispersion broad- 
ening time, tdm (ms), across one frequency channel, Av, is re- 
lated to the dispersion measure (DM) as 


^dm 


2 7tm e c 


(i 


v V 


DM 


8.3x10 6 ^^DM (54) 


v 3 

K MHz 


where m e is the mass of the electron, c is the speed of the light, 
and v \ , V 2 are the edges of the frequency channel. The dispersion 
measure, DM (pc cm -3 ), is obtained using the Cordes & Lazio 
(2001) NE2001 model. The same model provides the scattering 
measure, SM (kpc m -20/3 ), which allows to estimate the broad- 
ening time due to interstellar scattering as 


Tscat = 1000 


/ SM 

\292 


1.2 


i - 4.4 
^ V GHz 


(55) 


where d is the pulsar distance in kpc (Johnston et al., 1992; 
Stumer & Dermer, 1996). The last term of equation 53, T 2 ailDM , 
is an additional time broadening added when the sampling is 
performed for a DM value different from the real one. It corre- 
sponds to the fourth term of equation 2 in Dewey et al. (1985) 
and becomes important just for low period pulsars. 
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4. (°) 

4d O 

b s t O 

K d (°) 

dec st (°) 

dec e d (°) 

Gurv 

^Dewey 

Duty cycle 

Molonglo 2 

- 

- 

- 

- 

-85.0 

20.0 

0.62 

0.4 

0.03 

Green Bank 2 

- 

- 

- 

- 

-18.0 

90.0 

0.32 

0.7 

0.03 

Green Bank 3 

15.0 

230 

-15 

15 

- 

- 

0.41 

0.75 

0.03 

Parkes2 

- 

- 

- 

- 

-90.0 

0.00 

0.90 

0.75 

0.03 

Arecibo 2 

40 

66 

-10 

10 

9.50 

25.0 

0.54 

1.0 

0.05 

Arecibo 3 

38 

66 

-8.1 

8.2 

5.00 

26.5 

0.66 

0.7 

0.05 

Parkes 1 

-92 

20 

-4 

4 

- 

- 

0.41 

0.6 

0.03 

Jodrell Bank 2 

-5 

105 

-1.3 

1.3 

- 

- 

0.50 

0.8 

0.03 

Parkes MB 

-105 

52 

-6.03 

6.35 

- 

- 

0.98 

0.9 

0.05 

Swinburne 

-100 

50 

4.5 

30 

- 

- 

0.87 

1.0 

0.05 


Table 3. Estimated survey parameters. Respectively from the left to the right column are indicated: longitude start & end, latitude start & end, 
declination start & end, new survey efficiency, Dewey scaling factor, and pulsar duty cycle, defined as the pulse width over the period and used in 
the computation of the intrinsic pulse width W (Equation 53). 


The sky temperature at frequencies other than 408 MHz is 
obtained as: 

/ 408 MHz\ 2 ' 6 

Tsky(VMHz) = Tsky.408 • (56) 

\ VmHz / 

Tables 2 and 3 list all the radio telescope and detector charac- 
teristics of the surveys we took into account. Some of the survey 
parameters in the literature listed as average values have been 
re-evaluated using the above mentioned prescription. 

The scaling of the radiometer equation was motivated by the 
uncertainties related to the Dewey formula, because of flux os- 
cillations due to scintillation. The scintillation is caused by the 
turbulent variation of the interstellar medium that the pulsar light 
has to cross before reaching the observer. The consequence is 
an oscillation (scintillation) of the pulsar flux that can introduce 
spurious detections of pulsars with a flux lower than the survey 
threshold or that can cause the non-detection of pulsars with a 
flux higher than the survey threshold. So we scaled the S threshold 
level in order to take into account possible spurious detections 
or missed detections due to scintillation. S threshold should not be 
lower than the flux of the weakest pulsar of the survey. A rea- 
sonable estimate is to employ the average of the low-flux tail of 
the pulsars of the survey. 

In the ATNF database we can count how many pulsars fall 
within a survey boundary, how many would match the survey 
visibility criterion (flux > S threshold), and how many of these pul- 
sars have really been observed by the survey. The comparison of 
the ratios of the radio flux recorded for each pulsar to the min- 
imum visible flux S m [ n in its direction provides an estimate of 
the Dewey scaling factor. The scaling values ebewey are given in 
Table 3 and the distribution of the flux ratios is shown in Figures 
A.l, A.2, and A.3 (right plots) for each survey (only the ratios 
below 10 are displayed to focus near the visibility threshold). 

Then, for each survey, we derived the ratio between the num- 
ber of pulsars really detected by the survey and the total number 
of observable ATNF pulsars (the sum of the detected ones plus 
those that match the position and flux survey criteria but were 
not detected). We consider this last ratio as the new survey effi- 
ciency, 6 surv , i.e. the percentage of pulsars detected by the survey 
with respect to all the detectable ATNF pulsars in the survey re- 
gion. The new efficiency 6 surv is listed, for each Survey, in Table 
3. 

By using the newly estimated survey parameters listed in 
Tables 2 and 3, and by using the radiometer equation 52, the 
number of real pulsars that meet the visibility criteria of our sur- 
veys is 1430 (ATNF database, January 2012). We use this num- 
ber and the number of simulated pulsars that match the same cri- 


teria to scale the visible component of the simulated y-ray pulsar 
population in Equation 50. 


10. Results 

10.1. Detection statistics 

Table 4 indicates, for each model, the numbers of NSs that 
passed the radio and/or y visibility criteria and their compari- 
son with the LAT detections after 2 years of observations. The 
number of radio visible pulsars in the simulation has been scaled 
to the 1430 ATNF radio pulsars that passed the same selection 
criteria. The scale factor of 0.136 is required to match the sim- 



log(Per) (s) log(Per) (s) 


Fig. 6. Number density of the visible radio pulsars as a function 
of period and period derivative. The left and right plots respec- 
tively show the simulation and observed data with the same grey 
scale saturating at 25 star/bin and the same visibility criteria. The 
rising grey line marks the slot-gap death line. The declining grey 
lines mark the iso-magnetic lines at 10 7 , 10 8 , and 10 9 Tesla. 


ulated and observed radio samples and has been applied to all 
star counts quoted hereafter, in particular to the y-ray simu- 
lated samples. This scale factor implies a NS birth rate of ~ 3.7 
NS/century over the last 1 Gyr. The choice of radiative efficien- 
cies driven by a reasonable agreement in the S y D 2 evolutions 
with characteristic age shows that the wide beams produced in 
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PS R( Radio v y) 

PS -^(Radio A y) 

P S Ry only 

PS ^(Radio A y)/ PS Ry all 

PC 

1431 

3.6 

0.5 

0.87 

SG 

1508 

75 

79 

0.49 

OG 

1496 

123 

66 

0.65 

OPC 

1524 

107 

94 

0.53 

LAT 

\ 

25 

30 

0.45 


Table 4. For each model and for the observed dataset, we give from left to right the scaled numbers of pulsars visible in the radio or y-ray band, 
in the radio and y-ray bands, in y-rays only, and the fraction of radio loud objects in the y visible sample. All the data refer to two years of LAT 
observations. 


the intermediate-high (SG) and outer models provide enough de- 
tections to account the LAT findings. The low-luminosity narrow 
PC beam fails in predicting the LAT detection number and the 
fraction of radio-quiet objects because of the large overlap be- 
tween the y-ray and radio beams. 

10.2. Comparison of the total simulated and observed 
samples 

Figures 6 and 7 show the comparison between the simulated dis- 
tributions in the P - P diagram for the radio visible component 
and for the y-ray visible population for each model. The simu- 
lated distributions reasonably describe the observed samples and 
are in nice agreement with the same distributions obtained by 
Takataetal. (2011). The simulated radio population is able to 
describe the observed P - P distribution for the fastest rotators 
that are likely to sustain substantial y-ray emission and represent 
the LAT pulsar population. The PC model reproduces poorly 


1/5 



Fig. 7. Number density of the visible y-ray pulsars obtained for 
each model as a function of period and period derivative. The 
linear grey scale saturates at 1.5 star/bin. The pink triangles and 
green dots show the radio-loud and radio-quiet LAT pulsars, re- 
spectively. The rising grey line in the slot-gap subplot marks 
the slot-gap death line. The declining grey lines mark the iso- 
magnetic lines at 10 7 , 10 8 , and 10 9 T. 


the observed population. Both the SG and OG models over pre- 
dict the number of middle aged y-ray pulsars and under predict 
the number of young y-ray objects. Of those, the OG shows the 
poorer description of the data; the core of the distribution is too 
close to the pulsar death line and it lacks energetic pulsars. The 


OPC y- visible population best describes the observed P and P of 
the LAT population with a core centred on the observed objects 
and tails that cover the overall dispersion. 

Figure 8 compares the total simulated populations and its y- 
ray sub-sample to the observed total sample of radio and/or y-ray 
visible objects for key characteristics: period, period first time 
derivative, surface magnetic field, and distance. The simulated 
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sim. obs. Io 9( b ) CO 



0 5 10 15 

Distance (kpc) 


Fig. 8. Number distributions in period, period first time deriva- 
tive, surface magnetic field strength, and distance obtained for 
the whole populations of radio or y-ray visible pulsars in the 
simulations (light grey histogram) and in the LAT and radio sur- 
vey data (thick line). The ATNF radio sample has been restricted 
to the objects that pass the same position and flux selection crite- 
ria as in the simulation. The slot-gap model has been used as an 
example for the y-ray simulation. The dark shaded histograms 
show the distributions of the gamma active subsample of the 
whole simulated population. The abundances of simulated ob- 
jects at low P, high P, and high B are dominated by y-ray active 
pulsars. The excess of energetic and nearby simulated objects 
reflects the set of assumptions adopted for the birth distributions 
to provide a better match to the LAT data. 


spin period distribution is too broad to describe the observed 
proportion between the number of intermediate period objects 
(~ 50 ms) and the wings of the distribution. The range of spin 
periods is well covered and well centred, but we lack simulated 
objects in the 0.3 -1.0 second range. The simulated distributions 
in P, B , and D are all shifted to an excess of young, energetic, 
and nearby pulsars compared to the observed ones. This results 
from the choice of birth characteristics and NS intrinsic charac- 
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teristics (A/jvs, Rns, and I formulation) that emphasised nearby 
and high -E objects while preserving the bulk of the radio distri- 
butions. This choice has been made a posteriori to minimise the 
lack of high-E objects discussed in section 10.3. Nevertheless, 
the discrepancies observed in Figures 6 and 8 are not only due 
to the choice of birth distributions, but also to a radio model ill 
adapted to explain the observed radio population at the high- 
est E s. Whereas this would be problematic to study radio beam 
models, the reasonable representation at P < 500 ms and the ex- 
cess of objects with P > 3 x 10~ 15 s/s, B > 10 8 T, and D < 4 
kpc, where most of the LAT pulsars are found, supports the study 
of y-ray models. The necessity of an improved radio model is a 
result of this paper and its formulation, beyond the purpose of 
this study, will be the subject of future work. In the histograms 
shown in Figure 8, the total distributions are dominated by the ra- 
dio sample since the y-ray pulsars’ contribution is much smaller. 

10.3. The spin-down power 

Figures 9 and 10 compare the distributions in spin-down power 
and characteristic age for the LAT and the y- visible simulated 
pulsars. All models are significantly lacking simulated pulsars 
with spin-down power E > 3 x 10 28 W and characteristic age 
^char < 100 kyr. Additionally, all the models over-predict the 
number of low E pulsars and favour a pulsar population older 
than that observed by the LAT. The difference in shape of the 



26 28 30 32 26 28 30 32 
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Fig. 9. Spin-down power distributions obtained for each model 
for the visible y-ray pulsars. Pink and green refer to radio-loud 
and radio-quiet fractions of the total population, respectively. 
The LAT distribution (in black) has been scaled to the total 
number of visible pulsars for each model to ease the compari- 
son and show the relative lack of young energetic pulsars with 
E> 3 x 10 28 W. 


observed and simulated histograms suggests that the E inconsis- 
tency is not due to a simple scale mismatch, but to a deficiency 
in modelling the pulsar evolution: even by scaling the spin down 
power upward none of the models would be able to describe the 
observed distribution. 

Even though the PC model fails to produce enough visible 
gamma-ray pulsars because its narrow beam is under luminous 


and rarely visible, its evolution with E or age is less skewed to 
old age than the high-altitude SG or the outer-gap models. 

The OG model provides the poorest description of the y- 
ray evolution. A strong evolution with age is predicted by the 
classical formulation of the OG because the gap size is con- 
trolled by the amount of X-rays emitted by the stellar surface 
heated by the back-flow of primary charges returning from the 
gap (self-sustaining OG model). The strong evolution driven by 
this feedback is apparently not supported by the LAT data. The 
OPC model gives slightly better results but still fails to predict 
enough high E pulsars. The similarity of the E profiles obtained 
for the SG and OPC models shows that the relative lack of ener- 
getic y-ray pulsars is not related to the number of visible hemi- 
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Fig. 10. Age distributions obtained for each model for the visible 
y-ray objects. Pink and green refer to radio-loud and radio-quiet 
fractions of the total population, respectively. The LAT distri- 
bution (in black) has been scaled to the total number of visible 
objects for each model to ease the comparison and show the rel- 
ative lack of visible objects with age <100 kyr. 


spheres (two pole caustic SG, one pole caustic OPC), or to the 
evolution of the emission region with age within the open mag- 
netosphere (the emitting layer moves closer to the magnetic axis 
with increasing age in the OPC case while it remains near the 
last closed B line, but widens with age in the SG model). 

The under-prediction of high-i? visible y-ray pulsars is rather 
puzzling since they are the intrinsically brightest objects (high 
particle power and large fa) with the widest beams (large open 
magnetosphere and thin gaps emitting near the closed field lines) 
sweeping widely across the sky. The problem affects all the mod- 
els, so its origin does not depend much on the emission pat- 
tern or the luminosity trend with E. For instance, the luminos- 
ity evolution of the OPC model was constructed to agree with 
the LAT data, yet the deficit of energetic y-ray visible pulsars 
is still present. For a given luminosity the effective flux inter- 
cepted by the observer strongly depends on the gap thickness. 
For E > 10 28 W, the OPC gap width is 10 to 100 times smaller 
than the SG one, concentrating the photons in sharp caustics that 
remain visible to large distances and over many aspect angles, 
yet both the OPC and SG models produce a deficit of high E pul- 
sars in a rather similar way. The discrepancy is also insensitive to 
the relative orientation of the radio and y-ray beams since both 
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radio-loud and radio-quiet simulated pulsars are missing at high 
E. Nor is the problem related to the sensitivity horizon since all 
the models over-predict the fainter objects at low E. By testing 
different population configurations we have tried to understand 
which pulsar parameter has the largest impact on the high E tail 
of the y-ray sample. Different birth distributions in period and 
magnetic field have been tested. Decreasing the birth spin pe- 
riod in order to increase the fraction of very young and energetic 
pulsars (section 3.1) yields a very small gain in the number of 
y- visible energetic pulsars. Scanning the allowed domain of in- 
trinsic luminosities (e.g., SG A parameter) also failed to produce 
an increase in the young, energetic population. 


26 28 30 32 

RLRQ aLAT ^ 

Fig. 11. Spin-down power distributions obtained for each model 
for the visible y - ray objects for a different set of birth distri- 
butions: a Gaussian of width 200 ms centred at 0 for the peri- 
ods; the sum of two Gaussians in log 10 i? NS [Tesla], both 0.6 in 
width, respectively centred at 8.4 and 9.1, and with an amplitude 
ratio 1:7/12 for the surface magnetic fields; and the Paczynski 
(1990) surface density in the Galaxy. Pink and green refer to 
radio-loud and radio-quiet fractions of the total population, re- 
spectively. The LAT distribution, scaled to the total number of 
visible objects, is plotted as a black contour. 


A different choice of M N s, Rns , or / would shift the sim- 
ulated distributions horizontally in E, but would not alter their 
shape. The range of acceptable masses and radii given in 
Lattimer & Prakash (2007) limits an increase in the moment of 
inertia of the stars (thus E) to within a factor of 2 or 3 beyond 
our present choice. This is too small a factor to address the lack 
of high-E pulsars in the simulations. One of the tested config- 
urations, illustrated in Figure 11, shows how the lack of high- 
E pulsars remains, even after choosing a more energetic pulsar 
population at birth and a much broader distribution for their birth 
places across the Galaxy, Paczynski (1990), as used in previous 
work e.g., Gonthier et al. (2004) and Takata et al. (2011). 

Despite the stronger bias to energetic objects at birth adopted 
in Figure 11 as compared to Figure 9, the lack of high-E y - ray 
pulsars is less severe in the latter. This is due to the much larger 
fraction of births occurring in the inner Galaxy for the population 
shown in Figure 9. Because of the constraints on the supernova 
rate in the Galaxy, we cannot significantly increase the number 
of recent births, but the distribution provided by the HII region 
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Fig. 12. Number density of the visible radio and/or y-ray pul- 
sars in the Milky Way (polar view). The left and right plots re- 
spectively show the simulation and observed data with the same 
logarithmic gray scale saturating at 100 star/bin and the same 
visibility criteria. The cyan contour outlines the region where 
simulated SG y-ray pulsars are detectable. The cyan triangles 
show the location of the LAT pulsars. The yellow dot marks the 
Sun. 
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Fig. 13. Distance distributions obtained for each model for the 
visible y-ray objects. Pink and green refer to radio-loud and 
radio-quiet fractions of the total population, respectively. The 
LAT distribution (in black) has been scaled to the total number 
of visible objects for each model to ease the comparison and 
show the relative overabundance of nearby objects for the PC 
and SG models and under abundance of nearby objects for the 
OG and OPC models. 


profile concentrates a larger fraction of the recent births in the in- 
ner Galaxy, within the LAT visibility horizon. So, the E problem 
seems related both to the birth location and spin-down evolution 
of the pulsars. 

It is possible that the magnetic obliquity a decreases with 
age, as suggested by Young et al. (2010). First, the solid angle 
swept by the pulsar beam would decrease as a gradually de- 
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Fig. 14. Distribution of the y-ray beaming factors obtained for each model as a function of the spin-down power. Pink and green 
dots refer to radio-loud and radio-quiet y- ray visible objects, respectively. The solid and dotted lines give the best exponential fit to 
the radio-loud and radio-quiet data points, respectively. For the PC case, a unique fit was applied to the whole sample of radio-loud 
and radio-quiet objects. 


creases with age, so pulsars detected originally will later become 
invisible. Second, a has an impact on visibility through the gap 
width. This is illustrated by the difference in the E histograms 
obtained for the OG and OPC cases. They share the same emis- 
sion pattern, but the OG gap width depends on a while the OPC 
gap width is just proportional to E. Another speculative expla- 
nation would be a slower evolution of the dipole spin-down for 
very young and energetic objects, within the first 100 kyr. This 
hypothesis would need to be justified on the basis of theory. 


10.4. Spatial distribution in the Milky Way 

Figure 12 shows a polar view of the spatial density of visible 
radio or y-ray pulsars in the Galaxy, resulting from the birth lo- 
cation profile described in section 3.2. 

The majority of these pulsars are bom within the solar circle 
and the y-ray visibility contour agrees well with the Galactic re- 
gion where the bulk of the LAT pulsars have been detected. The 
radio visibility horizon is closer in the simulation than in reality, 
but the y-ray visibility horizon spans the right distance range. 
The visibility is therefore not the primary cause for the lack of 
high -E y-ray predictions discussed in the previous section. 

Figure 13 gives the distance distributions of the visible y-ray 
pulsars. An interesting trend, observed in all the models except 
for the PC, is that the radio loud to radio quiet ratio increases up 
to 4 or 5 kpc, and decreases down to zero at larger distance, im- 


plying that we lose the radio emission with distance more rapidly 
than we lose the y-ray signal. 

Since the measurement of the LAT pulsar distances is often 
affected by large uncertainties, the comparison between models 
and data in Figures 12 and 13 should be taken with care. 

10.5. The beaming factor fa 

The beaming factor (Watters et al., 2009) is defined as 

* = 4 (57) 

where L y , D , and {vF v ) are respectively the pulsar y-ray lumi- 
nosity, distance, and average energy flux. The beaming factor is 
the ratio of the total energy flux radiated over the 4n sr solid an- 
gle swept by the pulsar beam, after one complete rotation, to the 
phase-averaged energy flux observed at a given f 0 b s angle. Its 
value depends both on the intrinsic solid angle of the emission 
beam, on the beam inclination, and on the amount of energy that 
it contains. From Equations 57 and 46, the beaming fraction is 
calculated for each simulated light-curve as: 

_ f sin Cdl, jy n{<p, ()d<p 

Jo 2 n ' 

2 fo »(£>bs,0)<# 

Figure 14 shows the behaviour of the beaming factor as a 
function of E. For each model we have fitted the trend for the 
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PC 




Fig. 15. Distribution of the y- ray luminosities obtained for each model as a function of the spin-down power. Pink and green dots 
refer to radio-loud and radio- quiet y-ray visible objects, respectively. The LAT luminosities (black circles and crosses for radio-loud 
and radio-quiet objects respectively) have been derived using the energy-flux measurement and the fa value estimated from the fit 
to the simulated data for the particular spin-down power and radio-loud or quiet state of the LAT pulsars. The dotted line shows the 
100% efficiency boundary. 


radio-loud and radio-quiet components of the population by us- 
ing an exponential function. In the PC case there are too few vis- 
ible y-ray pulsars, both radio-loud or radio-quiet, to fit the evo- 
lution of the beaming factor with E for each type separately. The 
PC fa(E ) has been evaluated by merging the samples and fitting 
the global trend. The best-fit functions (/q,rl(^) and fa,RQ(E)) 
are indicated in Figure 14. 

For all models, the small decrease of fa with age (decreas- 
ing E) is due to the shrinking of the polar cap as the pulsar slows 
down. In the PC case, both the radio-quiet and radio-loud fa 
values are very dispersed and very small because of the colli- 
mated PC beam. For the SG, OG, and OPC, the fa distribution 
of the radio-quiet population component appears always more 
dispersed than the radio-loud one and it spans lower fa values, 
while the radio-loud objects exhibit higher and more highly con- 
strained beaming factors. Radio-quiet pulsars are generally seen 
at large | a-£ | impact angle where the y-ray caustics are fainter, 
so fa can reach low values. 

The SG case shows a minimal change in beaming factor with 
age for both the radio-loud and radio quiet pulsars because emis- 
sion from the bright caustic can be seen over most £ directions in 
the sky and because this model predicts a strong off-peak emis- 
sion. The SG radio-quiet and radio-loud beaming factors are cen- 
tred around fa = 1 . 

In the OG and OPC cases, we note a pronounced dispersion 
in fa , over 1 or 2 orders of magnitude, for the radio-quiet pul- 
sars. The OG and OPC models share the same emission pattern 


(phase-plot), thus the dispersion covers the same range of values. 
In the OPC case, the beaming factor increases up to E ~ 10 28 W 
and then stays constant around 0.8. Since all the OPC simulated 
pulsars at a given E have the same y-ray luminosity (by construc- 
tion), the observed spread in the fa values reflects the spread in 
beam flux as seen from different perspectives. It amounts typi- 
cally to less than a factor of 2 for radio-loud objects and more 
than one order of magnitude for radio-quiet ones. As the pulsars 
age in the OG and OPC models, there is an increasing separation 
of the gamma-ray and radio beams on the sky as the gamma-ray 
beam shrinks towards the the spin equator, producing a greater 
number of radio-quiet pulsars with small fa . 

In the outer gap models, the core of the fa distributions 
is consistent with the beaming factor obtained by Takata et al. 
(201 V) fa~ 0.4. 

10.6. y-ray luminosity trend with E 

Figure 15 shows, for each emission model, the evolution of the 
y-ray luminosity with the spin-down power and its comparison 
with the LAT results. The luminosity of the LAT pulsars has 
been computed from the measured pulsed flux using Equation 
57 with a beaming factor fa(E) obtained from the best fit plotted 
in Figure 14, according to their radio-loud or radio-quiet state. 

The observed evolution is roughly predicted by all the mod- 
els. Given the large dispersion in both the data and model predic- 
tions, the luminosity trend with E cannot be used to discriminate 
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between the gap models. In the OPC case, the L y (E) evolution is 
a built-in assumption of the model chosen to follow the observa- 
tions. 

The comparison with the LAT population indicates that the 
PC model is not luminous enough to account for the observed 
pulsars. Because of the large radiative efficiency (increased 
power in the gap), the SG luminosity reasonably follows the 
LAT data and the SG population best describes the observed 
trend. Since the y-ray emission is sustained by the particles gen- 
erating the same polar cap electromagnetic cascade, L y in both 
the PC and SG models follows the same trend, steepening from 
L y oc E l/2 to L y cc E with decreasing E when the pulsar puts out 
most of its spin-down power into y-rays. This trend is predicted 
but not yet observed because of the large dispersion in the LAT 
data points and large uncertainties in LAT pulsar distances. It is 
possible, however, that a more pronounced break to lower y-ray 
luminosities is required at low E. 

The OG luminosity evolution shows a different behaviour 
for high and low spin-down values. At E < 10 28 W, the gap 
width quickly saturates to a constant value which covers about 
three quarters of the open field volume. Only objects with large 
obliquities a remain visible and their luminosity scales linearly 
with E. They exhibit a small dispersion that is not consistent with 
the LAT data. At higher E values, the LAT pulsars fall within the 
range of predicted luminosities but they exhibit less dispersion 
than predicted by the model. 

The L y oc E 0 - 5 evolution predicted by the OG and SG models 
at high E is driven primarily by the evolution of the Goldreich- 
Julian current across the open magnetosphere. This is true if the 
feedback between particle acceleration and electrical screening 
from the cascading yields a rather stable maximum energy for 
the primaries and a stable fraction of this energy is radiated away 
in the cascade. A steeper evolution ( L y oc E) is expected for both 
polar-cap and outer-gap accelerators when the electrical screen- 
ing becomes inefficient and the gap fills a large part of the open 
magnetosphere. The LAT data in Figure 1 5 suggest a stronger lu- 
minosity evolution at E > 10 29 W than proposed by the current 
models. This conclusion appears to be robust because it applies 
to the very different radiation patterns tested in the OG and SG 
cases and because the beaming factors that have been used to de- 
rive the LAT luminosities show little evolution and little scatter 
from one pulsar to the next for E > 10 29 W for both models. 

One can note a larger scatter in the SG and OG luminosities 
plotted in Figure 15 than in the beaming factors for the same E 
range. This is due to the intrinsic variation in gap width result- 
ing from the variety of NS properties, amplified by the fact that 
the luminosity is proportional to the gap width cubed. The scat- 
ter in the SG luminosity distribution is driven by the spread in 
period and stellar magnetic field at each value of E. The scatter 
in the OG luminosity further builds on a strong dependence of 
the gap width with obliquity. In that sense, looking forward to 
a time when more precise distance estimates are obtained and 
when tighter constraints on the gap location in the outer mag- 
netosphere (from phase-resolved spectroscopy and light-curve 
studies) provide more reliable beaming factors for each pulsar, 
the dispersion in the luminosity plot can teach us about the diver- 
sity of young neutron stars that compose the LAT pulsar sample. 

10.7. Fractions of y-ray loud and radio-loud pulsars 

Figure 16 illustrates the change in visibility distance for pul- 
sars of different ages (E) and emission types (loud or quiet in 
the radio and y-ray bands). A similar figure has been shown in 


Watters & Romani (2011) to study the evolution of the pulsar 
visibility horizon. Data points for the LAT pulsars have been 
overlaid, but one should note that very few distance estimates 
exist for the new LAT pulsars that have been found through blind 
periodicity searches, because of the lack of radio dispersion mea- 
sures. Given the large uncertainties in distance, the agreement 
between the SG and OPC predictions and the LAT data is rea- 
sonable. 

Figure 17 shows the fraction of radio loud pulsars in the cu- 
mulative distribution of y-ray visible pulsars with E larger than 
the plotted value, and Figure 1 8 shows the fraction of y-loud pul- 
sars in the cumulative distribution of radio-visible pulsars above 
the given E (see Ravi et al. 2010). 

They jointly illustrate the high probability of detecting both 
the radio and y-ray beams in LAT objects with E > 10 30 W, 
in contrast to the predictions of all the models. This fraction 
remains high at all E values for the PC prediction, at variance 
with the data, because the radio and y-ray beams are produced 
in nearby regions of the magnetosphere. 

Figure 16 shows that the LAT visibility horizon for radio- 
loud and radio-quiet objects shrinks similarly with pulsar age 
(decreasing E ), so we lose both types at the same rate. The radio- 
loud to radio-quiet ratio evolves little for E < 10 30 W, so the 
radio-loud fraction in the whole y-ray sample (Figure 1 7) flattens 
with age to about 0.5 in good agreement with the LAT data. 

The value of this fraction is controlled by the effective differ- 
ence in flux sensitivity for the detection of a radio-loud or radio- 
quiet pulsar in the LAT data. Detecting a y-ray pulsation after 
phase-folding with the radio ephemerides requires fewer photon 
counts and trials than for blind periodicity searches. 




mod. obs. log(Edot) (W) log(Edot) (W) 


Fig. 17. Spin-down power evolution of the fraction of radio loud 
pulsars in the cumulative distribution of y-ray visible pulsars 
with E larger than the plotted value. The thin lines give the sim- 
ulation results for each model. The thick line gives the fraction 
evolution in the LAT sample. 


The use of the first pulsar catalogue sensitivity map for the 
radio-loud objects (Abdo et al., 2010) and of the blind-search 
sensitivity map for the radio-quiet objects (Dormody et al., 
201 1), both scaled to 2 years, brings an excellent agreement be- 
tween the model predictions and data at E < 10 3 ° W. The use 
of the catalogue map (presenting the lowest flux thresholds) for 
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Fig. 16. For each model, the linear gray scale shows the number density of the radio- or y- visible population, with a saturation at 
10 star/bin. Contours are given for the radio-only visible (cyan), y-only visible (green), and radio plus y visible (red) objects at 1% 
of the maximum number density for each type. The red triangles and green dots mark the location of the radio-quiet and radio-loud 
LAT pulsars, respectively. 


both types of pulsars lowers the predicted fraction by a factor 
of 2. Figure 16 shows that the y-ray visibility horizon shrinks 
much more rapidly than the radio one with pulsar power. So y- 
ray pulsars become undetectable faster than radio pulsars and 
the y-loud fraction in the radio-visible sample falls continuously 
down to a few percent at low E in Figure 18. 

Figure 17 shows that the radio-loud fraction in the LAT sam- 
ple becomes 1 at high E while the SG and OPC predictions re- 
main flat or decrease. The outer-gap predictions decline at high 
E because of a relative increase in the number of radio-quiet y- 
ray pulsars seen at large distances (see Figure 16 for the OPC). 
We checked that the SG prediction behaves similarly if we in- 
crease its radiative efficiency to match that of the OPC model 
at high E , so the radio-loud deficit is rather independent of the 
1-pole versus 2-pole radiation pattern from the outer regions. 
The open magnetosphere widens with E and wide beams are 
produced near the edge of the open volume by the efficiently 
screened, thin gaps. One would thus expect a larger overlap be- 
tween the radio and y-ray beams (i.e. an increase in the radio- 
loud to radio-quiet ratio) as E increases. Another effect must 
overcome this trend. It is related to the detectability of the large 
reservoir of faint radio-quiet pulsars in the outer magnetosphere 
gap models. The latter exist at large | a - impact angles and 
they dominate in potential number. The y-ray caustic emission 
extends to larger | a - angles as E increases, so the number of 
potentially visible y-ray pulsars (radio-quiet) increases. Then, as 
the dimmer parts of the caustic emission intercepted at high im- 


pact angle gradually passes the sensitivity threshold as the lumi- 
nosity increases, a larger fraction of the potential reservoir be- 
comes y-ray visible. In other words, the radio-loud probability 
increases because of the widening of the radio and y-ray beams, 
but the flux detectability of the large reservoir of faint radio-quiet 
objects increases even more and the net result is that the radio- 
loud to radio-quiet ratio can remain constant or decrease at high 
E. We checked, by lowering the luminosity of the SG and OPC 
models or by increasing the flux sensitivity threshold, that the 
radio beam widening effect becomes dominant when it is harder 
to detect faint radio-quiet objects. 

The LAT visibility is good enough for the simulations to pre- 
dict a small fraction of radio-loud objects at high E , at variance 
with the LAT data. The use of higher flux thresholds for y-ray 
detection would alleviate this deficit, but it would significantly 
deteriorate all the other observable distributions. The predicted 
deficit is robust against different gap locations and extents in the 
outer magnetosphere (2 pole or 1 pole emission, infinitely thin 
emitting layer for OPC and emission across a gap thinning with 
increasing E for the slot gap, emission above the null surface 
or reaching to lower altitudes). It is also robust against the gap 
width estimation (which impacts the caustic extent) since the SG 
and OPC gap width values we obtained differ by 1 0 to 30 at high 
E. The pronounced discrepancy suggests that the assumed radio 
beam is too narrow at young ages. A broader beam of radio emis- 
sion at higher altitude has also been suggested by Manchester 
(2005) and Ravi et al. (2010). 
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mod. obs. log(Edot) (W) log(Edot) (W) 

Fig. 18. Spin-down power evolution of the fraction of y-loud pul- 
sars in the cumulative distribution of radio-visible pulsars with 
E larger than the plotted value. The thin lines give the simulation 
results for each model. The thick line gives the fraction evolution 
in the LAT sample. 


We have tested flat distributions in a and f instead of 
isotropic ones. The flat distributions give more weight to small 
angles, therefore bias the samples toward more numerous radio- 
quiet pulsars. Adopting isotropic distributions in a and £ results 
in lower fractions of both radio and y-ray, with even lower pre- 
dictions at medium and high E for all models compared to the 
predictions shown in Figures 17 and 18. A rapid decrease of the 
magnetic obliquity, over a timescale of 1 Myr, has been sug- 
gested by Young et al. (2010) from their study of radio pulse 
widths. The discrepancy we find here applies to ages less than 
30 kyr, so magnetic alignment is unlikely to play a key role in 
reconciling the observed overabundance of young radio and y- 
ray pulsars and the model predictions. 

11. Summary 

The exceptional results obtained with the Fermi LAT telescope 
in the last few years offer the unique and exciting opportunity 
to constrain the physics of the pulsed y-ray emission by study- 
ing the early evolution of the pulsar population and its collective 
properties. We compared simulation predictions with Fermi LAT 
observations for this young ordinary pulsar population. 

We synthesised a radio and y-ray pulsar sample, assuming 
a core and cone model for the radio emission and y-ray emis- 
sion according to four gap models, the Polar cap (PC), Slot Gap 
(SG), Outer Gap (OG), and an alternative outer gap, the One 
Pole Caustic Model (OPC), that uses the OG beam geometry 
and a simple luminosity evolution with E consistent with the 
LAT data (Watters et al., 2009). We compared model expecta- 
tions and LAT data by applying y-ray and radio visibility crite- 
ria to our sample and by scaling it to the number of radio pulsars 
observed in the Milky Way. 

We found that the narrow beam of the low-altitude polar cap 
emission contributes at most a handful of pulsars in the LAT 
sample. The modelled luminosity is also too faint by an order of 
magnitude to account for the LAT data if one applies the average 
PC beaming factors we found for the given spin-down powers 
of the LAT pulsars. The large dispersion found in PC beaming 


factors, however, can substantially solve the luminosity discrep- 
ancy. We find that all the LAT pulsars are much more luminous 
(by 1 order of magnitude) than the PC expectations. Yet, there is 
a huge dispersion (1 or 2 orders of magnitude) in fa for the PC 
beams, so applying the average fa(E) trend to he LAT pulsars 
could be off by more than one order of magnitude in reality, so 
all the LAT points could go up and down by more than 1 decade 
in Figure 15 without problems. The wide beams from the outer 
gaps and slot gap models can easily account for the Fermi LAT 
detection number in 2 years, provided an increase of a factor of 
~ 10 of the standard slot-gap luminosity. The required increase 
may result from an enhanced accelerating electric field in the 
context of offset polar caps (Harding & Muslimov, 2011). The 
evolution of the enhanced SG luminosity with spin-down power 
is compatible with the large dispersion seen in the LAT data. 

We took into account the difference in the LAT flux sensi- 
tivity for detecting pulsed emission from radio-selected pulsars 
and for blind periodicity searches. The use of the two different 
sensitivity maps explained the almost equal amounts of radio- 
loud and radio-quiet pulsars found by the LAT. For all models, 
we found that the y-ray visibility horizon extends to compara- 
ble distances for radio-loud and radio-quiet pulsars as a func- 
tion of E , from 6 to 8 kpc at the highest powers down to 2 kpc 
for the least energetic LAT pulsars. The radio visibility horizon 
compares well with the y-ray horizon at high E , but it extends 
to much larger distances for less energetic pulsars, except for the 
rapidly evolving OG case for which the pulsars with E < 3 x 1 0 27 
W put 58% of their spin-down power into y-rays and remain vis- 
ible to 5 kpc. 

All the y-ray models fail to reproduce the high probability 
of detecting both the radio and y-ray beams at high E. The OPC 
prediction for the fraction of y-loud pulsars among the radio pul- 
sars is consistent with the radio and LAT data, but the model 
significantly under-predicts the fraction of y-ray pulsars that 
are radio-loud. The SG model also over-predicts the number of 
radio-quiet y-ray pulsars. These discrepancies may indicate that 
pulsar radio beams are larger than those we have modeled, ei- 
ther because they are intrinsically wider or because the emission 
occurs at higher altitude (Manchester, 2005), or both. The same 
conclusion has been argued by Karastergiou & Johnston (2007), 
that postulates emission over a wide range of emission heights 
rather than over a wide range of beam longitudes, and more re- 
cently by Ravi et al. (2010) and Watters & Romani (201 1) in the 
light of the Fermi observations. 

The beaming factor fa hardly evolves with E in the SG case. 
It is well constrained around 1 for both radio-loud and radio- 
quiet pulsars. In the OPC case, fa ~ 0.8 for radio-loud objects 
with E > 10 28 W, and it decreases by a factor of 2 for the less 
energetic objects detected by the LAT. In the OG case, fa de- 
creases from 1 to 0.3 with E decreasing down to 10 28 W and the 
evolution flattens around ~0.2 for lower powers. In all the mod- 
els, the beaming factor of radio-quiet pulsars follow the average 
trend found for radio-loud pulsars, but with a large dispersion 
than spans 1 or 2 orders of magnitude. 

The classical outer-gap model (OG) fails to explain many 
of the most important pulsar population characteristics, such as 
spin-down power distribution and luminosity evolution, whereas 
the outer-gap alternative (OPC), which is based on a simple 
scaling of the gap width with E~ l/2 , provides the best agree- 
ment between model predictions and data, as concluded by 
Watters & Romani (2011). This agreement relies on the very 
narrow gaps assumed in the OPC case. They are 10 to 100 times 
thinner than the values obtained for the SG for the same spin- 
down power, so the y-ray luminosity is concentrated in thin and 
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wide beams along the edge of the open magnetosphere. The OG 
model predicts a stronger luminosity evolution because it uses 
the polar cap heating by the returning particles to close the gap. 
The stronger evolution driven by this feedback is apparently not 
supported by the LAT data. Takata et al. (201 1) studied the evo- 
lution of the two layer OG luminosity as a function of E. Its 
result is consistent with the one we plot in Figure 15 for the OG 
model. The less pronounced dispersion observed in Takata et al. 
(201 1) is due to the choice of a fa - 1 for all the pulsars. 

All models studied here significantly under-predict the num- 
ber of visible y-ray pulsars seen at high E. This inconsistency 
does not depend on the modelling of the y-ray and radio visibil- 
ity thresholds. The discrepancy with the observations is signif- 
icant despite our choice of birth distributions skewed to young 
energetic pulsars, at slight variance with the constraints imposed 
by the total radio and y-ray pulsar sample observed. The fact that 
the four models have different y-ray luminosity evolutions and 
different beam patterns suggests a different cause for the discrep- 
ancy. Concentrating the birth location in the inner Galaxy less- 
ened but did not resolve the discrepancy. Further increasing the 
number of energetic pulsars near the Sun would conflict with the 
observed pulsar distances. The estimate of the visibility thresh- 
old in radio or y-ray flux is not at stake since all models over- 
predict the number of older, fainter, visible objects. 

The set of present results suggests that the observations re- 
quire rather luminous albeit thin gaps in the magnetospheres 
of young pulsars. It will be a challenge for models to match 
this behaviour. The impact of a magnetic alignment with age 
(Young et al., 2010), of an azimuthal variation of the accelerat- 
ing field, or the choice of different braking indices for the pul- 
sar spin-down may be important and will be included in future 
population studies to explore the origin of the scarcity of young 
energetic y-ray pulsars in the model predictions. 
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Appendix A: Radio survey sensitivity plots 
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Fig. A.l. Definition of the radio visibility criteria for the surveys Molonglo 2, Green Bank 2, Green Bank 3, and Parkes 2. For each 
survey we show the pulsar selection, classification, and counting in the redefined survey coordinates region (table 3), and the ratio 
between the pulsar fluxes and the threshold ones re-evaluated by taking into account the fudge factor defined in table 3. The object 
categories listed in the legend are: ATNF (black cross): pulsars inside the new survey coordinates region defined in Table 3 that are 
listed in the ATNF catalogue; ATNF vis (blue dot): pulsars that are visible according to the new survey parameters listed in Tables 2 
and 3 and/or listed in the ATNF catalogue; ATNF NaN (green dot): pulsars listed in the ATNF catalogue for which it is not possitti 
to evaluate the threshold flux; Detected (red circle): pulsars that have been detected according to the new survey parameters listed 
in Tables 2 and 3. 
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Fig. A.2. Definition of the radio visibility criteria for the surveys Arecibo 2, Arecibo3, Parkes 1, and Jodrell Bank 2. For each survey 
we show the pulsar selection, classification, and counting in the redefined survey sky region (table 3), and the ratio between the 
pulsar fluxes and the threshold ones re-evaluated by taking into account the fudge factor defined in table 3. The object categories 
listed in the legend are: ATNF (black cross): pulsars inside the new survey coordinates region defined in Table 3 that are listed in 
the ATNF catalogue; ATNF vis (blue dot): pulsars that are visible according to the new survey parameters listed in Tables 2 and 3 
Ifid/or listed in the ATNF catalogue; Detected(red circle): pulsars that have been detected according to the new survey parameters 
listed in Tables 2 and 3. 
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Fig. A.3. Definition of the radio visibility criteria for the surveys Parkes Multibeam and Swinburne. It is shown the pulsar selection, 
classification, and counting in the redefined survey sky region (table 3), and the ratio between the pulsar fluxes and the threshold 
ones re-evaluated by taking into account the fudge factor defined in table 3. The object categories listed in the legend are: ATNF 
(black cross): pulsars inside the new survey coordinates region defined in Table 3 that are listed in the ATNF catalogue; ATNF 
vis (blue dot): pulsars that are visible according to the new survey parameters listed in Tables 2 and 3 and/or listed in the ATNF 
catalogue; ATNF NaN (green dot): pulsars listed in the ATNF catalogue for which it is not possible to evaluate the threshold flux; 
Detected (red circle): pulsars that have been detected according to the new survey parameters listed in Tables 2 and 3. 
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